%WEATHER/PLOT Weather class plotting function.
%   w = PLOT(my_weather,time,altitude,type,fig) plots the data loaded into the
%   weather object MY_WEATHER at the specified time and altitude using figure
%   FIG.
%
%   TYPE can be one of the following:
%     'rain'          rain data
%     'rain il'       rain data interpolated between timesteps using linear
%                     interpolation.
%     'rain wind'     rain data and wind vectors
%     'rain wind i'   rain data and wind vectors interpolated between timesteps
%                     and datapoints using the INTERPN function.
%     'rain wind il'  rain data and wind vectors interpolated between timesteps
%                     using linear interpolation.
%     'rain wind il3' rain data and wind vectors interpolated between timesteps
%                     using linear interpolation, shown in 3d.
%     'hail'          hail data
%     'z wind'        vertical winds
%
%   elstonj 02-20-08
function plot(self,time,alt,type,fig)
%*=+--+=#=+--      UAS Severe Weather Simulation Softeware        --+=#=+--+=#*%
%          Copyright (C) 2013 Regents of the University of Colorado.           %
%                             All Rights Reserved.                             %
%                                                                              %
%    This program is free software: you can redistribute it and/or modify      %
%    it under the terms of the GNU General Public License version 2 as         %
%    published by the Free Software Foundation.                                %
%                                                                              %
%    This program is distributed in the hope that it will be useful,           %
%    but WITHOUT ANY WARRANTY; without even the implied warranty of            %
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the             %
%    GNU General Public License for more details.                              %
%                                                                              %
%    You should have received a copy of the GNU General Public License         %
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.     %
%                                                                              %
%           Jack Elston                                                        %
%           elstonj@colorado.edu                                               %
%                                                                              %
%*=+--+=#=+--                 --+=#=+--+=#=+--                    --+=#=+--+=#*%

settings = self.settings;
data = self.simdata;

if alt < self.simdata.z0
  %fprintf(1,'??? Input argument "altitude" must be >= %0.2f.\n\n',self.simdata.z0);
  %return;
  alt = self.simdata.z0;
end
if alt > self.simdata.zf
  %fprintf(1,'??? Input argument "altitude" must be <= %0.2f.\n\n',self.simdata.zf);
  %return;
  alt = self.simdata.zf;
end

%k=floor((alt-data.z0)/data.zs) +1;
k=min(find(data.z >= alt)); % work around to be able use non-uniform grids
t=floor((time-data.t0)./data.ts) +1;

x = data.x;
y = data.y;
z = data.z;

f = (time-data.t0) * data.frame.u;
g = (time-data.t0) * data.frame.v;

switch lower(type)

  case 'rain'

    for i = 1:size(t,2)
      figure(fig); hold off;

      if t(i) > data.tnum
        rain = data.rain(:,:,k,data.tnum);
      elseif t(i) < 1
        rain = data.rain(:,:,k,1);
      else
        rain = data.rain(:,:,k,t(i));
      end

      [cs,ch]=contourf(x+f(i),y+g(i),rain,0:0.001:0.009,'EdgeColor','none');
      set(gca,'color',[0,0,0.5625]);
      caxis([0 0.006]);

      title(strcat('Rain [kg/kg] time:',num2str((t(i)-1)*data.ts+data.t0),'s',' alt:',num2str(alt),'m'));
      xlabel('X Pos [m]');
      ylabel('Y Pos [m]');
      axis([data.x0 data.xf+data.frame.u*(data.tf-data.t0) data.y0 data.yf+data.frame.v*(data.tf-data.t0)]);
    end
		
  case 'z wind'

    for i = 1:size(t,2)
      figure(fig); hold off;

      if t(i) > data.tnum
        w = data.w(:,:,k,data.tnum);
      elseif t(i) < 1
        w = data.w(:,:,k,1);
      else
        w = data.w(:,:,k,t(i));
      end

      [cs,ch]=contourf(x+f(i),y+g(i),w,'EdgeColor','none');
      %set(gca,'color',[0,0,0.5625]);
      %caxis([-6 6]);

      title(strcat('Vertial Wind Speed [m/s] time:',num2str((t(i)-1)*data.ts+data.t0),'s',' alt:',num2str(alt),'m'));
      xlabel('X Pos [m]');
      ylabel('Y Pos [m]');
      axis([data.x0 data.xf+data.frame.u*(data.tf-data.t0) data.y0 data.yf+data.frame.v*(data.tf-data.t0)]);
    end
		
  case 'z wind cross'

		y_location = min(find(data.y >= 19000));

    for i = 1:size(t,2)
      figure(fig); hold off;

      if t(i) > data.tnum
					u = data.u(y_location,:,:,data.tnum) + data.frame.u;
					v = data.v(y_location,:,:,data.tnum) + data.frame.v;
					w = data.w(y_location,:,:,data.tnum);
					rain = data.rain(y_location,:,:,data.tnum);
      elseif t(i) < 1           
					u = data.u(y_location,:,:,1) + data.frame.u;
					v = data.v(y_location,:,:,1) + data.frame.v;
					w = data.w(y_location,:,:,1);
					rain = data.rain(y_location,:,:,1);
      else                        
					u = data.u(y_location,:,:,t(i)) + data.frame.u;
					v = data.v(y_location,:,:,t(i)) + data.frame.v;
					w = data.w(y_location,:,:,t(i));
					rain = data.rain(y_location,:,:,t(i));
      end

			%u = reshape(u,length(z),length(x));
			%w = reshape(w,length(z),length(x));
			u = reshape(u,length(x),length(z))';
			w = reshape(w,length(x),length(z))';
			rain = reshape(rain,length(x),length(z))';

			step = 8;

      %[cs,ch]=contourf(x+f(i),z,rain,0:0.001:0.009,'EdgeColor','none');
      [cs,ch]=contourf(x+f(i),z,rain,'EdgeColor','none');
      set(gca,'color',[1,1,1]);
      %caxis([0 0.006]);

			map = [1:-1/100:0]';
			colors = [map map map];
			colormap(colors);
			colorbar;

			hold on;

			h1=quiver((x(1:step:end)+f(i))./1,z,u(:,1:step:end),w(:,1:step:end));
      set(h1, 'Color',[0 0 0]);
      set(h1, 'AutoScaleFactor',0.5);
      set(h1, 'LineWidth',2);

      title(strcat('Vertial Wind Speed [m/s] time:',num2str((t(i)-1)*data.ts+data.t0),'s',' Y dist:',num2str(19000),'m'));
      xlabel('X Pos [m]');
      ylabel('Z Pos [m]');

			hold on;
			h = plot([x(1)+f(i) x(end)+f(i)],[300 300],'r','LineWidth',2);
			legend(h,['Transect']);

			axis([x(1)+f(i) x(end)+f(i) z(1)  z(end)+200]);
    end


  case 'z wind il'
	k = k+1;

    delta_t = mod(time,data.ts);
    for i = 1:size(t,2)
      figure(fig); hold off;

			spacing = 4;

      if t(i) > data.tnum - 1
        u = data.u(1:spacing:data.xnum,1:spacing:data.ynum,k,data.tnum) + data.frame.u;
        v = data.v(1:spacing:data.xnum,1:spacing:data.ynum,k,data.tnum) + data.frame.v;
        w = data.w(:,:,k,data.tnum);
      elseif t(i) < 1
        u = data.u(1:spacing:data.xnum,1:spacing:data.ynum,k,1) + data.frame.u;
        v = data.v(1:spacing:data.xnum,1:spacing:data.ynum,k,1) + data.frame.v;
        w = data.w(:,:,k,1);
      else
        u = (data.u(:,:,k,t(i)+1) - data.u(:,:,k,t(i)))/data.ts*(delta_t(i))+data.u(:,:,k,t(i)) + data.frame.u;
        u = u(1:spacing:data.xnum,1:spacing:data.ynum);
        v = (data.v(:,:,k,t(i)+1) - data.v(:,:,k,t(i)))/data.ts*(delta_t(i))+data.v(:,:,k,t(i)) + data.frame.v;
        v = v(1:spacing:data.xnum,1:spacing:data.ynum);
        w = (data.w(:,:,k,t(i)+1) - data.w(:,:,k,t(i)))/data.ts*(delta_t(i))+data.w(:,:,k,t(i));
      end

      [cs,ch]=contourf(x+f(i),y+g(i),w,'EdgeColor','none');
      set(gca,'color',[128,255,128]./255);
      caxis([-15 15]);

      hold on;

      h1=quiver(x(1:spacing:data.xnum)+f(i),y(1:spacing:data.ynum)+g(i),u,v);

      set(h1, 'Color',[0 0 0]);
      %set(h1, 'LineWidth',2);

      h1 = title(strcat('Wind [m/s] and Rain [kg/kg] time:',num2str((t(i)-1)*data.ts+data.t0),'s',' alt:',num2str(alt),'m'));
      set(h1, 'FontSize',16);

      title(strcat('Vertial Wind Speed [m/s] time:',num2str((t(i)-1)*data.ts+data.t0),'s',' alt:',num2str(alt),'m'));
      xlabel('X Pos [m]');
      ylabel('Y Pos [m]');
      axis([data.x0 data.xf+data.frame.u*(data.tf-data.t0) data.y0 data.yf+data.frame.v*(data.tf-data.t0)]);
    end

  case 'z wind journal'

    delta_t = mod(time,data.ts);
    for i = 1:size(t,2)
      figure(fig); hold off;

			step = 2;

      if t(i) > data.tnum - 1
        u = data.u(1:step:data.xnum,1:step:data.ynum,k,data.tnum) + data.frame.u;
        v = data.v(1:step:data.xnum,1:step:data.ynum,k,data.tnum) + data.frame.v;
        w = data.w(:,:,k,data.tnum);
      elseif t(i) < 1
        u = data.u(1:step:data.xnum,1:step:data.ynum,k,1) + data.frame.u;
        v = data.v(1:step:data.xnum,1:step:data.ynum,k,1) + data.frame.v;
        w = data.w(:,:,k,1);
      else
        u = data.u(:,:,k,t(i)) + data.frame.u;
        u = u(1:step:data.xnum,1:step:data.ynum);
        v = data.v(:,:,k,t(i)) + data.frame.v;
        v = v(1:step:data.xnum,1:step:data.ynum);
        w = data.w(:,:,k,t(i));
      end

      [cs,ch]=contourf(x+f(i),y+g(i),w,'EdgeColor','none');
      set(gca,'color',[1,1,1]);
      caxis([-15 15]);

			map = [0:1/100:1]';
			colors = [map map map];
			colormap(colors);
			colorbar;

      hold on;

      h1=quiver(x(1:step:data.xnum)+f(i),y(1:step:data.ynum)+g(i),u,v);

      set(h1, 'Color',[0 0 0]);
      set(h1, 'LineWidth',2);

      h1 = title(strcat('Wind [m/s] and Rain [kg/kg] time:',num2str((t(i)-1)*data.ts+data.t0),'s',' alt:',num2str(alt),'m'));
      set(h1, 'FontSize',16);

      title(strcat('Vertial Wind Speed [m/s] time:',num2str((t(i)-1)*data.ts+data.t0),'s',' alt:',num2str(alt),'m'));
      xlabel('X Pos [m]');
      ylabel('Y Pos [m]');
    end

  case 'rain il'
      
    delta_t = mod(time,data.ts);
    for i = 1:size(t,2)
      figure(fig); hold off;

      if t(i) > data.tnum - 1
        rain = data.rain(:,:,k,data.tnum);
      elseif t(i) < 1
        rain = data.rain(:,:,k,1);
      else
        rain = (data.rain(:,:,k,t(i)+1) - data.rain(:,:,k,t(i)))/data.ts*(delta_t(i))+data.rain(:,:,k,t(i));
      end

      [cs,ch]=contourf(x+f(i),y+g(i),rain,0:0.001:0.009,'EdgeColor','none');
      set(gca,'color',[0,0,0.5625]);
      caxis([0 0.006]);

      title(strcat('Wind [m/s] and Rain [kg/kg] time:',num2str((t(i)-1)*data.ts+data.t0),'s',' alt:',num2str(alt),'m'));
      xlabel('X Pos [m]');
      ylabel('Y Pos [m]');
      axis([data.x0 data.xf+data.frame.u*(data.tf-data.t0) data.y0 data.yf+data.frame.v*(data.tf-data.t0)]);
    end

  case 'wind'

    for i = 1:size(t,4)
      figure(fig); hold off;

			if (data.xnum > 4) && (data.ynum  > 4)

				if t(i) > data.tnum
					u = data.u(1:4:data.xnum,1:4:data.ynum,k,data.tnum) + data.frame.u;
					v = data.v(1:4:data.xnum,1:4:data.ynum,k,data.tnum) + data.frame.v;
				elseif t(i) < 1
					u = data.u(1:4:data.xnum,1:4:data.ynum,k,1) + data.frame.u;
					v = data.v(1:4:data.xnum,1:4:data.ynum,k,1) + data.frame.v;
				else
					u = data.u(1:4:data.xnum,1:4:data.ynum,k,t(i)) + data.frame.u;
					v = data.v(1:4:data.xnum,1:4:data.ynum,k,t(i)) + data.frame.v;
				end

				%h1=quiver(x(1:2:data.xnum)+f(i),y(1:2:data.ynum)+g(i),u,v);
				h1=quiver(y(1:4:data.ynum)+g(i),x(1:4:data.xnum)+f(i),u,v);

			else

				if t(i) > data.tnum
					u = data.u(:,:,k,data.tnum) + data.frame.u;
					v = data.v(:,:,k,data.tnum) + data.frame.v;
				elseif t(i) < 1
					u = data.u(:,:,k,1) + data.frame.u;
					v = data.v(:,:,k,1) + data.frame.v;
				else
					u = data.u(:,:,k,t(i)) + data.frame.u;
					v = data.v(:,:,k,t(i)) + data.frame.v;
				end

				data

				h1=quiver(y+g(i),x+f(i),u,v);

			end

      set(h1, 'Color',[0 0 0]);

      title(strcat('Wind [m/s] time:',num2str((t(i)-1)*data.ts+data.t0),'s',' alt:',num2str(alt),'m'));
      xlabel('X Pos [m]');
      ylabel('Y Pos [m]');
			if (data.xf ~= data.x0) && (data.yf ~= data.y0)
				axis([data.x0 data.xf+data.frame.u*(data.tf-data.t0) data.y0 data.yf+data.frame.v*(data.tf-data.t0)]);
			end
      hold off;
    end

  case 'rain wind'

    for i = 1:size(t,2)
      figure(fig); hold off;

      if t(i) > data.tnum
        rain = data.rain(:,:,k,data.tnum);
        u = data.u(1:2:data.xnum,1:2:data.ynum,k,data.tnum) + data.frame.u;
        v = data.v(1:2:data.xnum,1:2:data.ynum,k,data.tnum) + data.frame.v;
      elseif t(i) < 1
        rain = data.rain(:,:,k,1);
        u = data.u(1:2:data.xnum,1:2:data.ynum,k,1) + data.frame.u;
        v = data.v(1:2:data.xnum,1:2:data.ynum,k,1) + data.frame.v;
      else
        rain = data.rain(:,:,k,t(i));
        u = data.u(1:2:data.xnum,1:2:data.ynum,k,t(i)) + data.frame.u;
        v = data.v(1:2:data.xnum,1:2:data.ynum,k,t(i)) + data.frame.v;
      end

      [cs,ch]=contourf(x+f(i),y+g(i),rain,0:0.001:0.009,'EdgeColor','none');
      set(gca,'color',[0,0,0.5625]);
      caxis([0 0.006]);

      hold on;

      h1=quiver(x(1:2:data.xnum)+f(i),y(1:2:data.ynum)+g(i),u,v);
      set(h1, 'Color',[1 1 1]);

      title(strcat('Wind [m/s] and Rain [kg/kg] time:',num2str((t(i)-1)*data.ts+data.t0),'s',' alt:',num2str(alt),'m'));
      xlabel('X Pos [m]');
      ylabel('Y Pos [m]');
      axis([data.x0 data.xf+data.frame.u*(data.tf-data.t0) data.y0 data.yf+data.frame.v*(data.tf-data.t0)]);
      hold off;
    end

  case 'rain wind journal'

    for i = 1:size(t,2)
      figure(fig); hold off;

			step = 4;

      if t(i) > data.tnum
        rain = data.rain(:,:,k,data.tnum);
        u = data.u(1:step:data.xnum,1:step:data.ynum,k,data.tnum) + data.frame.u;
        v = data.v(1:step:data.xnum,1:step:data.ynum,k,data.tnum) + data.frame.v;
      elseif t(i) < 1
        rain = data.rain(:,:,k,1);
        u = data.u(1:step:data.xnum,1:step:data.ynum,k,1) + data.frame.u;
        v = data.v(1:step:data.xnum,1:step:data.ynum,k,1) + data.frame.v;
      else
        rain = data.rain(:,:,k,t(i));
        u = data.u(1:step:data.xnum,1:step:data.ynum,k,t(i)) + data.frame.u;
        v = data.v(1:step:data.xnum,1:step:data.ynum,k,t(i)) + data.frame.v;
      end

      [cs,ch]=contourf(x+f(i),y+g(i),rain,0:0.001:0.009,'EdgeColor','none');
      set(gca,'color',[1,1,1]);
      caxis([0 0.006]);

			map = [1:-1/100:0]';
			colors = [map map map];
			colormap(colors);
			colorbar;

      hold on;

      h1=quiver(x(1:step:data.xnum)+f(i),y(1:step:data.ynum)+g(i),u,v);
      set(h1, 'Color',[0 0 0]);
      set(h1, 'LineWidth',2);

      title(strcat('Wind [m/s] and Rain [kg/kg] time:',num2str((t(i)-1)*data.ts+data.t0),'s',' alt:',num2str(alt),'m'));
      xlabel('X Pos [m]');
      ylabel('Y Pos [m]');
      hold off;
    end

  case 'rain wind 3'
    x = self.simdata.x;
    y = self.simdata.y;
    z = self.simdata.z;
    t = self.simdata.time;

    u = interpn(x,y,z,t,self.simdata.u,x,y,z,time);
    v = interpn(x,y,z,t,self.simdata.v,x,y,z,time);
    w = interpn(x,y,z,t,self.simdata.w,x,y,z,time);
    rain = interpn(x,y,z,t,self.simdata.rain,x,y,z,time);
    hail = interpn(x,y,z,t,self.simdata.hail,x,y,z,time);

    figure(fig); hold off;
    f = (time-data.t0) * data.frame.u;
    g = (time-data.t0) * data.frame.v;

    [sx sy sz] = meshgrid(x,y,z);
    verts = stream3(x,y,z,u,v,w,sx,sy,sz);
    return
    sl = streamline(verts);
    iverts = interpstreamspeed(x,y,z,u,v,w,verts,.025);
    axis tight; view(30,30); daspect([1 1 .125])
    camproj perspective; camva(8)
    set(gca,'DrawMode','fast')
    box on
    streamparticles(iverts,35,'animate',10,'ParticleAlignment','on')
    %[cs,ch]=contourf(x+f,y+g,rain,0:0.001:0.009,'EdgeColor','none');

  case 'rain wind i'

    x = self.simdata.x;
    y = self.simdata.y;
    z = self.simdata.z;
    t = self.simdata.time;

    u = interpn(x,y,z,t,self.simdata.u,x,y,alt,time);
    v = interpn(x,y,z,t,self.simdata.v,x,y,alt,time);
    rain = interpn(x,y,z,t,self.simdata.rain,x,y,alt,time);
    hail = interpn(x,y,z,t,self.simdata.hail,x,y,alt,time);

    figure(fig); hold off;
    f = (time-data.t0) * data.frame.u;
    g = (time-data.t0) * data.frame.v;

    [cs,ch]=contourf(x+f,y+g,rain,0:0.001:0.009,'EdgeColor','none');
    set(gca,'color',[0,0,0.5625]);
    caxis([0 0.006]);
    axis([data.x0 data.xf+data.frame.u*(data.tf-data.t0) data.y0 data.yf+data.frame.v*(data.tf-data.t0)]);
    hold on;
    %colorbar;

    h1=quiver(x(1:2:data.xnum)+f,y(1:2:data.ynum)+g,u(1:2:data.xnum,1:2:data.ynum),v(1:2:data.xnum,1:2:data.ynum));
    set(h1, 'Color',[1 1 1]);
    title(strcat('Rain [kg/kg] time:',num2str(time),'s'));
    xlabel('X Pos [m]');
    ylabel('Y Pos [m]');


  case 'rain wind il'
      
    delta_t = mod(time,data.ts);
    for i = 1:size(t,2)
      figure(fig); hold off;

      if t(i) > data.tnum - 1
        rain = data.rain(:,:,k,data.tnum);
        u = data.u(1:2:data.xnum,1:2:data.ynum,k,data.tnum) + data.frame.u;
        v = data.v(1:2:data.xnum,1:2:data.ynum,k,data.tnum) + data.frame.v;
      elseif t(i) < 1
        rain = data.rain(:,:,k,1);
        u = data.u(1:2:data.xnum,1:2:data.ynum,k,1) + data.frame.u;
        v = data.v(1:2:data.xnum,1:2:data.ynum,k,1) + data.frame.v;
      else
        u = (data.u(:,:,k,t(i)+1) - data.u(:,:,k,t(i)))/data.ts*(delta_t(i))+data.u(:,:,k,t(i)) + data.frame.u;
        u = u(1:2:data.xnum,1:2:data.ynum);
        v = (data.v(:,:,k,t(i)+1) - data.v(:,:,k,t(i)))/data.ts*(delta_t(i))+data.v(:,:,k,t(i)) + data.frame.v;
        v = v(1:2:data.xnum,1:2:data.ynum);
        rain = (data.rain(:,:,k,t(i)+1) - data.rain(:,:,k,t(i)))/data.ts*(delta_t(i))+data.rain(:,:,k,t(i));
      end

      [cs,ch]=contourf(x+f(i),y+g(i),rain,0:0.001:0.009,'EdgeColor','none');
      set(gca,'color',[0,0,0.5625]);
      caxis([0 0.006]);

      hold on;

      h1=quiver(x(1:2:data.xnum)+f(i),y(1:2:data.ynum)+g(i),u,v);
      %set(h1, 'Color',[1 1 1]);

      set(h1, 'Color',[0 0 0]);
      set(h1, 'LineWidth',2);

      h1 = title(strcat('Wind [m/s] and Rain [kg/kg] time:',num2str((t(i)-1)*data.ts+data.t0),'s',' alt:',num2str(alt),'m'));
      set(h1, 'FontSize',16);

      xlabel('X Pos [m]');
      ylabel('Y Pos [m]');
      axis([data.x0 data.xf+data.frame.u*(data.tf-data.t0) data.y0 data.yf+data.frame.v*(data.tf-data.t0)]);
      hold off;
    end

  case 'rain wind il3'
      
    x = self.simdata.x;
    y = self.simdata.y;
    z = self.simdata.z;

    u = self.simdata.u(:,:,:,t);
    v = self.simdata.v(:,:,:,t);
    w = self.simdata.w(:,:,:,t);
    rain = self.simdata.rain(:,:,:,t);

    figure(fig); hold off;
    f = (time-data.t0) * data.frame.u;
    g = (time-data.t0) * data.frame.v;

    %slice(x, y, z, rain,[],[],125:1000:2375);
		%hold on;
		[X,Y,Z] = meshgrid(x,y,z);
		h1=quiver3(X,Y,Z,u,v,w);
		%hold off;

    title(strcat('Rain [kg/kg] time:',num2str(time),'s'));
    xlabel('X Pos [m]');
    ylabel('Y Pos [m]');

  case 'hail'

    for i = 1:size(t,2)
      figure(fig); hold off;

      if t(i) > data.tnum
        hail = data.hail(:,:,k,data.tnum);
      elseif t(i) < 1
        hail = data.hail(:,:,k,1);
      else
        hail = data.hail(:,:,k,t(i));
      end

      [cs,ch]=contourf(x+f(i),y+g(i),hail,0:0.001:0.009,'EdgeColor','none');
      set(gca,'color',[0,0,0.5625]);
      caxis([0 0.006]);

      title(strcat('Hail [kg/kg] time:',num2str((t(i)-1)*data.ts+data.t0),'s',' alt:',num2str(alt),'m'));
      xlabel('X Pos [m]');
      ylabel('Y Pos [m]');
      axis([data.x0 data.xf+data.frame.u*(data.tf-data.t0) data.y0 data.yf+data.frame.v*(data.tf-data.t0)]);
    end

  case 'particle'
  clf;
    x = self.simdata.x;
    y = self.simdata.y;
    z = self.simdata.z;
    t = self.simdata.time;

    numpoints = 50; 

    xi = self.simdata.x0:(self.simdata.xf-self.simdata.x0)/numpoints:self.simdata.xf;
    yi = self.simdata.y0:(self.simdata.yf-self.simdata.y0)/numpoints:self.simdata.yf;
    zi = self.simdata.z;

    u = interpn(x,y,z,t,self.simdata.u,xi,yi,zi,time);
    v = interpn(x,y,z,t,self.simdata.v,xi,yi,zi,time);
    w = interpn(x,y,z,t,self.simdata.w,xi,yi,zi,time);
    rain = interpn(x,y,z,t,self.simdata.rain,xi,yi,zi,time);
    velocity = sqrt(u.*u + v.*v + w.*w);

    %[sx sy sz] = meshgrid(xi,yi,self.simdata.z0);
    [sx sy sz] = meshgrid(self.simdata.xf,yi,1250);
    verts = stream3(xi,yi,zi,u,v,w,sx,sy,sz);
    sl = streamline(verts);
    view(3);

    iverts = interpstreamspeed(xi,yi,zi,u,v,w,verts,1);
    set(gca,'DrawMode','fast')
    box on
    streamparticles(iverts,5,'animate',5,'ParticleAlignment','on')

  case 'radar image'
    for i = 1:size(t,2)
      %figure(fig); hold off;

      if t(i) > data.tnum
        rain = data.rain(:,:,k,data.tnum);
      elseif t(i) < 1
        rain = data.rain(:,:,k,1);
      else
        rain = data.rain(:,:,k,t(i));
      end

		maximum = max(max(rain));
		gain = 9/maximum;
		rain = flipud(rain) .* gain;
		rain = uint8(rain);
		map = [  0   0   0;
		         0 255   0;
		        51 153   0;
					 255 255   0;
					 204 153  51;
					 255 153   0;
					 255  51   0;
					 204  51   0;
					 153  51   0;
		       255  51 255]./ 255;

		imwrite(rain, map, fig,'gif','TransparentColor', 0, 'WriteMode', 'overwrite');

		x_length = m_to_lon(self.simdata.xf-self.simdata.x0,self.simdata.offset.lat);
		y_length = m_to_lat(self.simdata.yf-self.simdata.y0);
		width = self.simdata.xnum;
		height = floor(self.simdata.ynum*y_length/x_length);

		command = sprintf('convert -adaptive-resize %dx%d\\! %s %s',width,height,fig,fig);
		system(command);
    end

  case default

    fprintf(1,'??? Plot type %s not defined.\n\n',type);

end
